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1.   INTRODUCTION 


In  1981-1983  a  computer  program  was  developed  which  solves  numerically 
two-dimensional  Euler  equations  using  the  Godonov  method  (tlef.  1).   The 
program  was  intended  primarily  for  turbomachinery  application,  but  with 
modification  to  the  boundary  conditions,  application  to  any  problem  which  is 
modeled  by  the  Euler  equations  is  possible. 

This  report  covers  one  year  of  extensive  testing  of  the  code  for 
different  flow  regimes.   This  testing  was  done  with  three  main  objectives: 

a)  To  gain  experience  in  simulation  of  the  va.rious  problems  related  to 
the  research  activity  in  the  Turbopropulsion  Laboratory, 

Naval  Postgraduate  School,  Monterey,  California. 

b)  To  test  code  accuracy  by  comparing  numerical  with  analytical 
solutions . 

c)  To  find  the  source  of  errors  for  various  numerical  simulations  and  th? 
basic  Godunov  method  in  order  to  reduce  these  errors. 

Some  of  the  results  of  these  experiments  with  the  code  were  reported  in 
other  papers  and  appear  as  appendices  to  this  report.  Other  cases  were 
simulated,  but  not  published,  and  will  be  reported  here.  Some  unsuccessful 
attempts  to  improve  the  code  accuracy  are  also  described. 

The  topics  studied  and  the  progress  made  are  described  in  the  following 
sections . 

1)  The  problem  of  the  gradual  opening  of  wave  rotor  passage. 

2)  Numerical  modeling  of  the  nonsteady  thrust  produced  by  intermittent 
pressure  rise  in  a  diverging  channel. 

3)  Numerical  techniques  for  wave  rotor  cycle  analysis  (A.  Mathur). 

4)  Development  of  a  grid  generation  program  in  order  to  generate  a  grid 
over  wedge-arc  configuration. 

Test  case  of  the  supersonic  flow  in  a  channel  with  arc,  wedge, 
wedge-arc,  sinsoidal  arc-bumps. 

6)  Modification  of  the  Godunov  method  to  include  shock  fitting  using  the 
oblique  shock  wave  theory. 

7)  Development  of  the  one-dimensional  code  based  on  Verhoff's  method. 

3)   Modification  of  the  program  for  cylindrical  geometry  problems,  and 
modeling  of  the  flow  in  CDTD. 

9)   Modeling  of  the  flow  in  a  detonation  engine. 

10)   Development  of  the  SELCO  Method. 

LI)   Grid  generation  for  CDTD  inlet 


2.    PROBLEMS  CONSIDERED 


2. 1   The  Gradual  Opening  in  Wave  Rotor  Passage. 

Influence  on  the  flow  pattern  of  the  gradual  passage  opening  in  the  wave 
rotor  was  studied  on  the  basis  of  a  numerical  simulation.   It  was  found  that, 
in  most  cases,  a  significant  volume  of  the  passage  will  have  rotational  flow, 
which  should  lead  to  the  mixing  between  the  driver  and  driven  gases.   In  some 
cases,  losses  will  occur  as  a  result  of  the  multiple  reflections  of  the  shock 
and  pressure  waves  from  the  passage  walls.  .  It  was  shown  that  the  interface 
between  driven  and  driver  gas  will  be  oblique  to  the  passage  walls  when  the 
passage  opens  gradually,  and  that  the  interface  can -retain  its  obliqueness  to 
the  walls.   The  results  for  the  rectangular  and  skewed  passages  are  reported 
in  Appendix  A. 

In  the  case  of  the  skewed  passage,  the  perfect  compression  of  the  driven 
gas  with  very  small  mixing  losses  could  be  achieved  if  the  channel  opens 
according  to  a  certain  formula  of  the  optimum  opening  velocity.   It  was  found 
that  the  velocity  of  opening  could  be  predicted  by  the  equation: 

vopen  =  -. T  (1) 

sin  Lsk 

where:   ^open  =  the  velocity  of  the  gradual  opening 

vsh   =  the  shock  wave  velocity  in  the  driven  gas 
^■sk   =  tne  angle  of  the  skewed  passage 
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Numerical  Modeling  of  the  Nonsteady  Thrust  Produced  by  Intermittent 
Pressure  Rise  in  a  Diverging  Channel. 


The  dynamics  of  the  expansion  of  detonation  products  in  a  diverging 
channel  were  investigated  numerically.   The  influence  of  particular  features 
in  the  expansion  process,  such  as  the  presence  of  reversed  flow  and 
propagation  of  hammer  shocks  on  the  production  of  thrust  were  examined. 
Sequential  expansions  of  detonation  products  were  also  modeled  and  It  was 
concluded  that  in  order  to  maintain  a  high  frequency  periodic  mode  of 
operation  for  propulsion  applications,  the  channel  should  be  refilled  with 
ambient  air  after  each  expansion.   The  influence  of  the  ratio  of  ambient  air 
to  detonation  products  volume  on  the  dynamics  of  the  thrust  production  and  on 
the  impulse  generated  during  the  expansion  are  also  reported.   The 
presentation  of  the  results  of  this  investigation  is  given  in  Appendix  B. 

2.3   Numerical  Techniques  for  Wave  Rotor  Cycle  Analysis. 

With  the  author's  participation,  A.  Mathur  developed  a  one-dimensional 
code  based  on  Godonov's  method  for  wave  rotor  cycle  inalysis.   Because  the 
rasolution  of  the  code  on  the  interfaces  was  poor,  it  was  decided  to  develop  a 
one-dimensional  code  based  on  the  random  choice  method.   Non-standard  boundary 


conditions  for  this  problem  were  developed.   A  paper  was  written  and  accepted 
for  publication  in  the  proceedings  of  the  AS ME  Symposium  on  Nonsteady  FLows. 
The  paper  is  attached  as  Appendix  C. 

2.4  Development  of  the  Grid  Generation  Program  for  Wedge-Arc  Configurations. 

The  grid  generation  routine,  which  was  used  for  the  cascade  and  wave 
rotor  simulations  did  not  allow  a  full  control  of  the  grid  parameters.   This 
control  was  especially  needed  for  the  simulation  of  the  test  problems,  since 
we  wanted  to  insure  that  the  grids  have  the  same  general  parameters  for  the 
different  wedge-arc  configurations  considered. 

A  program  which  allows  the  grid  generation  with  full  control  over  grid 
parameters  was  developed.   The  program,  with  an  example  of  the  grid,  is  given 
in  Appendix  D  was  delivered  to  the  Naval  Postgraduate  School  and  stored  in  an 
archive  file. 

2 . 5  Test  Case  of  the  Supersonic  Flow  in  a  Channel  with  Arc,  Wedge ,  Wedge-Arc, 
Sinusoidal  Arc-Bumps. 

Performance  of  the  Godunov  code  was  tested  for  numerous  supersonic  and 
subsonic  internal  flow  problems.   The  tasks  included: 

a)  An  attempt  to  determine  what  causes  the  nonsymetrical  behavior  for 
subsonic  flow  case  and  whether  it  is  possible  to  increase  the 
accuracy  of  the  basic  Godunov  code  by  changing  boundary  conditions. 

b)  A  comparison  of  the  numerical  solution  with  analytical  predictions 
for  supersonic  flow. 

For  subsonic  flows  it  .«7as  determined  that  the  main  source  of  errors  was  slop 
discontinuity  at  the  corners  of  the  "bump".  When  the  circular  arc  bump  was 
replaced  by  a  sinusoidal  arc  "bump"' ,  the  symmetry  of  the  subsonic  solution 
improved  substantially.   This  resulted  from  the  fact  that  in  the  case  of  the 
sinusoidal  arc  only  the  second  derivatives  are  discontinuous  at  the 
surface.  It  was  found  that  the  boundary  conditions  at  the  surface,  which  are 
used  by  the  basic  Godunov  method  (solution  of  the  Riemann  problem  between  the 
physical  point  and  it  mirror  image  with  the  same  pressure  and  density  but  with 
the  negative  velocity),  produce  the  most  consistent  results  for  various 
problems.   Results  for  this  simulation  are  shown  in  Figure  1.   The  grid  is 
shown  in  Figure  2. 

Efforts  to  improve  the  subsonic  flow  solution  were  discontinued  vhen  it 
-/as  realized  that  a  more  objective  criteria  than  symmetry  of  solution  was 
needed.   For  this  reason,  work  on  supersonic  flows  was  begun,  since,  for  these 
conditions,  some  analytical  solutions  are  available.   First,  the  method  for 
supersonic  flow  (M  =  2.2)  in  a  channel  with  a  10%  thick  circular  arc  "bump  was 
tried.  The  grid  for  this  case  is  shown  in  Figure  3.   In  Figure  4  we  show 
results  for  this  geometry  when  the  inlet  Mach  number  is  2.0.   Results  are 
shown  in  the  form  of  the  distributions  of  the  pressure  coefficient  over  the 
surface  of  the  lower  wall  of  the  channel.   Numerical  results  correspond  very 
well  to  the  analytical  in  this  case.   It  was  found,  however,  that  when  other 
flow  parameters  are  compared  with  the  analytical  solution,  there  was 
substantial  disagreement  between  the  resits.   Figure  5  shows  the  comparison  of 
the  pressure  coefficient,  entropy,  density,  velocity  and  tach  number  that  ver^ 


calculated  numerically  with  the  analytical  solution.   The  comparison  is  done 
for  supersonic  flow  with  M=1.6  over  5%  thick,  wedge-arc  "bump".   As  shown  in 
Figure  5,  the  errors  are  produced  mostly  on  the  shock  waves  and  that  the 
largest  is  the  relative  error  in  entropy.   The  fact  that  only  pressure  is 
simulated  accurately  was  confirmed  for  the  numerous  test  cases  in  which  the 
arc  thickness  and  the  wedge  angle  varied  in  the  wide  range  of  values. 

2.6  Modification  of  the  Godunov  Method  based  on  the  Shock  Fitting  using  the 
Oblique  Shock  Wave  Theory. 

A  shock  fitting  was  attempted  to  reduce  the  errors  on  the  oblique  shocks. 
The  procedure  for  fitting  was  as  follows: 

1)  Calculate  the  values  using  Godunov' s  method. 

2)  Using  oblique  shock  wave  theory,  calculate  the  exact  value  of  the 
parameters.  This  calculation  is  based  on  the  assumption  that  the 
pressure  is  computed  correctly  by  the  Godunov  method. 

3)  Substituting  values  behind  the  shock  for  the  exact  values.  It  was 
found  that  the  procedure  did  not  work  if  it  was  applied  as 
described  above.   Many  modifications  of  this  prodedure  were  tried. 
Some  involved  interpolation  of  the  values  ahead  and  behind  the 
shock.   The  last  approach  did  improve  the  accuracy  of  the  entropy 
calculation  but  introduced  oscillations. 

2.7  Development  of  the  One-Dimensional  Code  based  on  the  Verhoff  Met 'nod. 

A  new  formulation  of  the  Euler  equations  and  numerical  method  to  solve 
them  were  proposed  by  Verhoff  (Ref.  2).   A  basic  code  implementing  his  method 
for  one -dimensional  flow  was  developed.   In  Figure  6,  results  are  shown  for 
the  case  of  two  colliding  shock  waves  in  a  tube.   These  results  are  a 
reconstruction  of  the  results  presented  in  Ref.  1  in  Figure  5.   Comparison  of 
these  two  figures  shows  that  our  code  implements  the  basic  Verhoff  method 
correctly.   The  listing  of  the  program  is  presented  in  Appendix  E. 

We  were  unable  to  implement  the  shock  fitting  in  our  code.   Comparison 
between  the  analytical  solution  and  numerical  was  10-15%  in  error  for 
velocity,  pressure  and  density  in  the  region  of  the  shock  wave,  if  the  shock 
was  not  fitted. 

2.8  Modification  of  the  Program  for  Cylindrical  Geometry  Problems,  and 
Modeling  or  the  Flows  in  CDTD: " 

The  Godunov  code  was  modified  for  simulation  of  the  flow  with 
cylindrical  velocity.   The  modified  program  was  used  for  modeling  of  the  flow 
field  in  CDTD  device  (Ref.  3).  The  task  of  the  modeling  was  to  test  how  the 
rotational  component  of  velocity  will  influence  flow  field. 

In  Figure  7,  one  of  the  geometries  of  the  CDTD   is  shown  with  a 
computational  grid  covering  the  domain  of  integration.   In  Figure  8,  the 
Iso-Mach  lines  for  the  geometry  are  shown  for  the  inlet  flow  conditions 
indicated   in  figure  7. 


The  conclusions  of  preliminary  research  was  that  the  program  models  the 
flow  correctly  and  that,  for  the  tested  conditions , the  results  are  close  to 
those  given  by  the  program  developed  by  Hirsch  as  discussed  in  Ref.  4. 

2.9  Modeling  of  the  Flow  in  Detonation  Engine, 

The  flow  field  in  the  projected  combustion  chamber  of  a  detonation 
engine  was  simulated  numerically  using  the  Godunov  code.   The  task  of  this 
research  was  to  study  the  dynamics  of  the  discharge  of  the  high  pressure  gas 
after  detonation.   A  number  of  configurations  were  tried.   In  Figure  9,  an 
example  of  this  type  of  simulation  can  be  seen.   Pressure  in  the  detonation 
chamber  drops  from  15  atm  to  0.55  atm  in  0.00131  second.   It  was  important  to 
establish  that  time  since  it  is  one  of  the  limitations  on  the  detonation 
frequency.   In  this  example,  the  ability  of  the  program  to  simulate  large 
gradients  of  pressure  can  be  seen. 

2.10  Development  of  the  SELCO  Method. 

A  new  method  for  supersonic  flow  simulation  was  developed.   A  full 
description  of  the  method,  with  illustrations  of  its  application  to  the 
supersonic  flow  over  the  wedge,  is  given  in  Appendix  D.   This  method 
demonstrates  that  inaccurate  modeling  of  the  oblique  shock  waves  produced  by 
the  Godunov  method  is  the  result  of  the  obliqueness  of  the  shock  wave  with 
respect  to  the  edges  of  the  cells  of  the  computational  grid  covering  domain  of 
integration.   It  was  shown  that  only  when  the  shock  surface  is  parallel  to  the 
two  opposite  edges  of  the  cell  that  the  oblique  shock  can  be  accurately 
calculated. 

A  new  method  of  the  Local  Cell  Orientation,  -  SELCO,  was  proposed  to 
allow  local  reorientation  of  the  cells  in  the  vicinity  of  the  shock  waves. 
The  efficiency  of  the  SELCO  method  was  demonstrated  for  the  simulation  of  the 
oblique  shock  waves  in  the  supersonic  flow  over  the  wedge. 

2.11  Grid  Generation  for  CDTD  Inlet. 

For  a  given  distribution  of  the  pressure  measurement  ports  on  the  CDTD 
blade  surface,  a  computational  grid  was  developed.   Several  variants  of  the 
CDTD  geometry  were  implemented.   However  all  grids  which  were  generated  proved 
to  lead  to  large  computational  errors.  This  work  was  not  continued  due  to  time 
constraints  on  the  project. 

3.     Conclusions 

The  principal  findings  of  the  program  of  code  testing  wece  as  follows: 

1.  The  code  has  a  robust  ability  to  analyse  unsteady  time  dependent 
flows  and  provide  qualitative  information  on  complex  shock  patterns 
and  pressure  fields. 

2.  With  attention  to  boundary  conditions  and  particularly  grid 
structure  Lt  is  possible  to  reproduce  analytical  flow  solutions  with 
great  accuracy. 


3.  The  solutions  produced  are  sensitive  to  grid  orientation.   The  grid 
must  be  re-orientated  locally  along  shock  or  discontinuity  planes. 

4.  Care  must  be  used  in  selecting  the  grid  distribution  on  areas  of 
rapid  curvature  change. 

Despite  the  program  sensitivity  under  certain  conditions,  with  skilled 
use  it  can  be  most  instructive  in  establishing  wave  and  snock  structures 
in  complex  flow  situations.   A  significant  example  is  the  unsteady 
boundary  of  Section  2.1  where  in  it  was  possible  to  generate  a  minimum 
loss  passage  angle  for  a  wave  rotor  device. 
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The  Problem  of  Gradual  Opening  in  Wave  Rotor  Passages 


S.  Eidelman* 
Naval  Postgraduate  School,  Monterey,  California 


I  In-  influence  on  ihe  flow  pattern  of  the  gradual  passage  opening  in  Ihe  wave  roior  is  studied  on  the  basis  of 
numerical  simulation.  It  is  found  (hat,  in  most  cases,  significant  volume  of  ihe  passage  will  have  rotational  flow, 
which  should  lead  to  the  mixing  between  the  driver  and  driven  gases.  In  some  cases,  losses  will  occur  also  as  a 
result  of  Ihe  multiple  reflection  of  Ihe  shock  and  pressure  waves  from  ihe  passage  walls.  It  is  -shown  lhai  the 
interface  between  driven  and  driver  gas  will  be  oblique  to  the  passage  walls,  when  Ihe  passage  opens  gradually, 
and  the  interface  can  retain  its  obliqueness  lo  ihe  walls. 


Introduction 

WAVE  rotors  are  devices  that  use  wave  propagation 
through  the  fluid  in  a  rotor  passage  to  transfer  energy 
from  one  fluid  to  another'  or  to  transfer  energy  from  one 
fluid  to  rotor  shaft  and  another  fluid. :  Wave  rotors,  wave 
engines,  wave  pressure  exchangers,  wave  equalizers,  and 
Comprex1  are  all  based  on  the  same  idea  of  energy  exchange 
in  the  unsteady  waves,  and  ihe  topic  of  the  present  study  is 
relevant  to  all  oft  hese  devices. 

Wave  rotors  offer  the  potential  for  significant  im- 
provements in  air-breathing  engine  propulsion  cycles  because 
they  can  be  self-cooling,  since  both  high-pressure  gas  and  cold 
low-pressure  gas  use  the  same  passage  for  alternate  periods  of 
lime  in  the  cycle.  Combining  a  wave  rotor  with  conventional 
turbomachinery  components  shows  promise  of  significant 
reduction  in  specific  fuel  consumption  without  weight  or  size 
penalties. 

The  principles  of  operation  of  wave  rotor  devices  and  their 
commerical  applications  can  be  found  in  Refs.  1-3.  For 
completeness,  the  general  scheme  of  a  wave  pressure  ex- 
changer will  be  described,  as  illustrated  in  Fig.  1.  One  gas 
(driver)  at  high  pressure  is  used  to  compress  a  second  gas 
(driven).  The  process  is  arranged  to  occur  in  tube-like 
passages  with  irapezoidal  cross  section  located  on  the  per- 
iphery of  a  rotating  drum  or  rotor.  The  compression  is 
achieved  successively  in  each  rotor  passage  by  means  of 
compression  waves  or  shock  waves  generated  by  ihe  entering 
dri\cr  gas.  Ihe  compressed  driven  gas  is  drawn  off  from  the 
end  of  the  passage  when  it  aligns  with  an  outlet  port.  The 
driver  gas  then  undergoes  a  series  of  expansions  lo  a  lower 
pressure  and  is  scavenged  out  by  freshly  inducted  driven  gas  at 
approximately  the  same  pressure  level.  This  fresh  "charge"  is 
then  compressed  by  the  high-pressure  driver  gas  and  the  cycle 
repeats  itself.  .Steady  rotation  of  the  drum  sequences  the  ends 
of  the  passages  past  stationary  inlet  ports,  outlet  ports,  and 
endwalls.  This  establishes  unsteady  but  periodic  flow  pro- 
cesses within  the  rotating  passages  and  essentially  steady  flow 
in  the  inlei  and  outlet  ports.  The  passage  may  be  oriented 
axially  as  in  fig.  I  or  at  a  stagger  angle.  In  general,  wave 
machines  used  as  pure  pressure  exchangers  (e.g.,  "Com- 
prex"') have  usually  axially  oriented  passages,  while  those 
with  staggered  passages  may  drive  a  shaft,  since  shaft  work 
extraction  is  possible  with  this  laiter  configuration. 
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In  the  design  of  the  wave  rotor,  it  is  very  desirable  to 
determine  the  optimal  ratio  between  the  width  and  length  ol 
the  single  passage  of  the  rotor.  Usually  only  two  parameters 
are  evaluated  to  determine  this  ratio:  skin  friction  losses  and 
bypass  losses.  The  number  of  passages  on  the  rotor  should  be 
minimal  to  minimize  skin  friction  losses.  On  the  other  hand, 
the  passage  should  be  narrow  compared  to  the  port  widths  to 
reduce  flow  bypass  between  inlet  or  outlet  ports.  The  transient 
process  of  the  passage  opening  or  closing  (as  the  passage  end 
moves  across  a  port  or  moves  from  a  port  to  be  closed  b>  the 
endwall,  respectively)  usually  is  not  considered  in  the  design. 
It  is  generally  assumed  that  the  passage  opening  or  closure 
occurs  instantaneously. 

Pearsoir  tried  to  take  into  account  the  gradual  opening  ot 
the  passage,  assuming  that  the  air  in  the  passage  is  eompres  -ed 
in  a  series  (usually  three)  of  discrete  compression  waves  w  Inch 
converge  and  ultimately  merge  to  form  a  shock  wave.  This 
allowed  him  to  design  a  complicated  wave  machine  cycle  for  a 
rotor  using  relatively  short  passages.  However,  since  ihe 
technique  was  one-dimensional  it  could  not  reveal  ihe 
peculiarities  of  this  essentially  two-dimensional  flow  and 
would  be  valid  only  for  very  weak  waves. 

in  the  present  study,  by  means  of  numerical  modeling,  we 
will  examine  in  real-time  how  the  gradual  opening  influences 
the  wave  formation  in  the  wave  rotor  passages  and  how  thai 
should  affect  the  rotor  design. 

Model 

The  assumptions  involved  in  the  numerical  simulations  are 
described  as  follows. 

The  flow  in  each  passage  of  the  wave  rotor  is  unsteady  and 
periodic.  At  the  same  time,  the  flow  through  ihe  port-,  is 
(ideally)  sieady.'  :  The  peripheral  width  ol  the  port  is  usually 
equal  to  several  widths  ot  a  single  passage,  and  herein  n  is 
assumed  that  the  flow  in  ihe  inlet  or  outlet  port  remains 
stationary  when  the  passage-end  encounters  the  port.  For  ihis 
reason  the  region  of  the  port  is  not  included  in  the  com- 
putational domain  shown  in  lie.  2. 

The  time-dependent  process  ai  the  passage-end  translating 
across  the  region  of  the  inlet  or  outlet  port  will  be  referred  .is 
the  "gradual  opening"  of  the  passage.  The  passage  opening 
process  will  be  referred  as  "instantaneous,"  when  the  as- 
sumption is  made  that  the  passage  instantaneously  opens  io 
the  port  area  and  is  subjected  immediately  to  the  steady  flow 
conditions  at  ihe  port. 

1 1  is  assumed  that  the  flow  is  inviscid  and  can  be  modeled  by 
i he  Fuler  equations.  Ihe  unsteady  two-dimensional  Fuler 
equations  can  be  w  rillen  in  conservation  law  form  as 

d(J      dh'      dG 

—  +  -     +  —  =0 
dt       dX      d  Y 
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Here  p  is  the  density,  u  and  v  the  velocity  components  in  the  X 
and  Y  coordinate  directions,  p  the  pressure,  and  7  the  ratio  of 
specific  heats.  The  energy  per  unit  of  volume,  e,  is  defined  by 


i-pi 


u2  +  v 2  \ 
—J-) 


where  t  =  P/(y-l)p  is  the  internal  energy.  We  look  for  the 
solution  of  the  system  of  equations  represented  by  Eq.  (I)  in 
the  computational  domain  shown  in  Fig.  2  in  time  t  with  the 
following  conditions  at  the  domain  boundaries:  a)  solid  wall 
along  segments  1-3  and  2-4,  b)  outflow  along  segment  3-4, 
and  c)  inlet  along  segment  1-2. 

It  is  assumed  that  initially  at  time  /  =  0,  the  passage  of  (he 
wave  rotor  is  filled  with  stationary  gas  at  ambient  conditions. 
When  instantaneous  opening  of  the  wave  rotor  passage  was 
simulated,  it  was  assumed  that  at  time  1  =  0,  the  flow  at  the 
inlet  1-2  was  equal  to  the  steady  flow  in  the  port.  When 
gradual  opening  of  the  passage  was  simulated,  at  time  /  =  0, 
the  inlet  was  closed  and  solid  wall  boundary  conditions  were 
imposed  at  the  inlet  1-2.  Then,  this  boundary  condition  was 
gradually  replaced  by  the  flow  condition  at  the  inlet  port.  The 
length  uncovered  to  the  inlet  port  region,  where  the  solid  wall 
boundary  conditions  were  replaced  by  the  inlet  port  con- 
ditions, was  determined  using  the  elapsed  time  and  the 
velocity  of  the  passage  relative  to  the  inlet  port. 

The  Godunov  method  was  used  to  obtain  a  numerical 
solution  of  Eq.  (1)  with  the  described  boundary  and  initial 
conditions.  Details  of  the  implementation  of  the  method  and 
boundary  conditions  are  given  in  Ref.  4. 

The  flowfield  was  modeled  for  the  rectangular  passage  with 
a  width  of  0.02  m  and  a  length  of  0. 12  m.  The  grid  covering 
the  computational  domain  of  the  passage  is  shown  in  Fig.  2. 

Results  and  Discussion 

The  following  initial  conditions  were  assumed  for  the  air 
initially  in  the  passage: 

/>„  =  latm     p0=1.2kg/mJ      Uo  =  0     Vo  =  0 

The  driver  gas  entering  through  the  port  at  the  left-hand  end 
was  assumed  to  have  the  following  properties: 

P,^  =  1.8atm     p,,=  1.81  kg/m3      l/,,  =  150m/s      V,,=0 

These  conditions  correspond  to  a  practical  situation  in  a  wave 
rotor  when  a  passage  filled  with  a  quiescent  fresh  charge  of  air 
at  ambient  conditions  encounters  an  inlet  port  supplying  hot, 
high-pressure  driver  gas. 

It  (he  assumption  is  made  that  the  passage  instantaneously 
opens  and  is  subjected  to  the  conditions  of  the  inlet  flow  at  the 
left  boundary  1-2  (see  Fig.  2),  then  a  perfectly  one- 
dimensional  How  pattern  should  develop  in  the  passage.  The 
results  of  ihe  modeling  of  these  conditions  are  presented  in  t he 
form  of  pressure  contours  at  progressively  larger  time  sieps  111 
Fig.  3.  Time  /  =  ()  corresponds  to  the  moment  when  the 
passage  opens.  Instantaneous  opening  of  the  passage  is  seen  in 
Fig.  3  to  lead  to  an  immediate  formation  of  a  shock  wave  and 
subsequent  propagaiion  in  the  passage  from  ihe  left  to  the 
right.  The  flow  condiiions  at  the  inlet  port  are  matched  10  the 
parameters  of  the   rarefaction   wave   which   effectively   will 
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Fig.  I     Wave  rotor  operation  scheme. 
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Fig.  2     Computational  domain. 
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Fig.  3     F'volulion  in  time  of  Ihe  shock  wave  formed  after  instant; 
ons  passage  opening. 


cancel  the  rarefaction  wave  moving  toward  ihe  inlet  at  ihe 
end  of  the  passage.  For  t his  reason  the  flow  in  the  passage 
no  additional  discontinuity.  This  situation  is  typical  ffl 
wave  rotor,  where  flow  conditions  at  the  ports  are  chosei 
such  a  way  that  waves  do  not  propagate  from  the  passe 
into  inlet  or  outlet  ports.  The  flow  in  the  ports  will  iheret 
remain  steady.  In  the  case  shown  in  Tig.  3,  ihe  shock  \\ 
was  found  to  propagate  with  a  velocity,  P\,h=446  m/s, 
interface  with  a  velocity  f ',„  =  150  m/s,  and  all  ihe  paramc 
examined  were  confirmed  accurately  using  one-dimcnsn 
gasdynamics  relationships. 

Most  of  the  approaches  used  in  wave  rotor  design  are  ba 
on  the  assumption  lhal  ihe  waves  are  one-dimensiona 
nature.  When  the  gas  is  compressed  by  a  weak  shock  wi 
for  instance,  assumptions  are  made  that,  l)  the  proces 
isoentropic,  2)  the  hot  and  cold  gases  are  strictly  separatee 
a  planar  interface,  and  3)  the  How  is  everywhere  irrotatio 
This  leads  to  ihe  very  high  efficiencies  projected  tor 
compression.  If  the  passage  is  wide  enough  so  thai  vise 
effects  can  be  neglected,  ihis  model  ol  the  compression  in 
wave  rotor  passage  is  realistic,  bin  only  for  instaiiianc 
passage  opening  or  for  a  very  long  passage.  Results  preset 
below  illustrate  how  ihe  gradual  passage  opening  a  I  feels 
compression  process. 
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In  Fig.  4,  the  pressure  and  Mach  number  contours  in  the 
passage  are  shown  at  a  sequence  of  times  tor  the  case  of  a 
gradual  opening  to  the  inlet  port  with  a  velocity  of  100  m/s. 
Parameters  lor  the  gas  in  the  passage,  before  opening  begins, 
and  in  the  inlet  port  are  (he  same  as  lor  the  case  of  in- 
stantaneous opening.  The  dynamics  of  the  Mow  development 
seen  in  Fig.  4a  is  very  different  from  that  shown  in  Fig.  3. 
First,  curved  pressure  waves  radiate  from  the  mitial  small 
openine  appearing  at  the  lower  corner  of  the  inlet  on  the  left 
side  of  the  passage  (r  =  0.044  ms  in  Fig.  4).  Subsequently, 
these' waves  reflect  from  the  upper  wall  of  the  passage  and  at 
time  f  =  O.I25  ms  have  formed  a  row  of  compression  waves 
which  have  approximately  straight  fronts  normal  to  the  wall 
of  the  passage.  Initially,  compression  waves  of  this  kind 
occupy  a  small  part  of  the  region  with  disturbed  gas.  In  the 
region  well  behind  the  front  of  this  quasi-one-dimensional 
propagation,  compression  waves  are  curved  and  are  the  result 
of  the  interaction  between  the  waves  reflecting  back  and  forth 
between  the  lower  and  upper  walls  of  the  passage  and  the  new 
waves  created  by  the  progressive  opening  of  the  port  to  the 
passage.  Since  it  is  possible  to  see  in  Fig.  4  for  the  times 
/=0.084  and  0.166  ms,  respectively,  the  flow  behind  the 
quasi-one-dimensional  region  is  highly  rotational  and  is 
relatively  constant  pressure  in  the  X  direction.  The  passage 
became  fully  opened  only  at  /  =  0.2  ms.  At  this  time  com- 
pression waves  are  propagating  along  the  length  of  the 
passage  and  the  pressure  rises  gradually  from  1  atm  at  the 
right  end  to  1.8  atm  at  the  left  end.  The  front  of  converging 
compression  waves  will  eventually  become  a  shock  wave,  but 
this  will  occur  at  some  later  time  and  outside  the  com- 
putational domain. 

It  was  concluded  from  the  case  presented  in  Fig.  4  that  a 
passage  length-to-width  ratio  of  6  will  lead  to  very  high 
mixing  losses  and  nonuniform  and  inefficient  compression, 
since  in  this  case  the  region  of  rotational  flow  occupies  half  of 
the  passage  volume.  Additional  losses  will  be  produced 
because  of  the  highly  rotational  flow  within  each  of  the  two 
gases. 

In  Fig.  5,  pressure  and  velocity  contours  arc  shown  for  the 
case  ol  gradual  opening  of  the  passage  with  a  velocity  of  200 
m/s.  lull  opening  of  the  passage  in  this  case  occurred  at 
/  =  ().l  ins.  A  curved  shock  wave  is  formed  at  /  =  0.044  ms. 
i This  shock  wave  (see  Fig.  5a)  partially  reflects  from  the  upper 
wall  of  the  passage  and  then,  gradually  converging  with  its 
main  Ironl,  forms  an  almost  planar  shock  front  at  the  time 
f  =  0.252  ms.  However,  even  then  the  flow  is  highly  rotational 
behind  the  shock  front,  and  the  region  of  rotational  flow 
occupies  one-third  of  the  passage  volume.  In  this  region  the 
gas  velocity  increases  from  iV/  =  0.3  at  the  upper  wall  to 
kw-0.52  at  the  lower  wall. 

When  i he  velocity  ol  the  passage  opening  is  increased  to  500 

!  m/s  (sec  I  ig.  6)  the  passage  becomes  fully  open  at  /  =0.04  ms. 

Because  of  the  last  opening  ol  the  passage,  the  shock  wave  at 

jthe  tune  i  =0.04  ms  is  less  curved  and  only  a  small  fraction  o\ 

I  the  shock  front  reflects  from  the  upper  wall.  This  reflected 

I  part  oi  the  shock   trout  converges  with  the  main   trout   at 

/  =  I).I72   ms.    From   that   time  on,   the   (low  pattern  in   the 

^passage     is     mostly     one-dimensional     with     a     small     and 

(weakening  region  ot  rotational  (low  behind  the  main  front. 

Nevertheless,  even  lot  ilus  case,  with  passage  length/passage 

width     3,  there  will  be  very  high  mixing  losses,  with  a  large 

part  ol  i he  passage  volume  subjected  to  rotational  How. 

In  mder  to  siudv  how  the  strength  of  the  shock  wave  in- 
fluences i  he  flow  pal  tern  developing  in  the  passage,  an  ad- 
Jiiional  case  was  simulated  where  the  parameters  ot  the  driver 
ias  ai  the  inlet  port  area  were  increased;  namely,  to 

/',,  =  2.X5  atm     U,,  -  283  m/s     p(/  =  4kg/nv' 

\s  in  the  previous  ease,  the  parameters  ol  the  driver  gas  were 
■hosen  in  such  a  way  that  waves  do  nol  propagate  back  into 
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arc  shown  in  fig.  7.  The  velocity  of  the  passage  opening  was 
I  =250  m/s.  It  was  concluded  from  the  results  presented  in 
Fig.  7  that  an  increase  in  the  initial  shock  strength  leads  io 
stronger  reflections  and  substantial  increases  in  the  How 
rotation.  This,  m  turn,  will  lead  to  increased  losses  because  of 
mixing  of  the  driver  and  driven  gases.  The  pattern  of  multiple 
reflections  of  (he  shock  wave  from  the  upper  and  lower  walls 
ol  i he  passage  can  be  seen  in  I  ig.  7a.  At  time  I- 0.073  ms  the 
shock  wave  reflected  from  the  upper  wall  is  propagating 
toward  the  lower  wall.  Part  ot  the  reflected  shock  wave  trout 
converges  with  the  main  shock  wave,  and  part  begins  to 
reflect  from  the  lower  wall  (r  =  0.109  ms).  I  he  multiple 
reflection  weakens  the  reflected  shock,  but  the  reflection 
between  the  walls  of  the  passage  continues  and  can  be 
followed  tor  /  -0.211  and  0.245  ms.  It  is  clear  t hat  lor  these 
conditions  even  passages  with  a  length-to-width  ratio  ol  6  will 
have  substantial  mixing  losses  because  ot  the  rotational  I  low. 

A  dil  leicui  situation  is  tound  tor  the  passage  opening  to  the 
inlet  port  ol  ihe  dnvcn  gas.  At  ihis  port,  the  pressure  and 
velocity  ol  the  (driver)  gas  in  the  passage  should  be  matched 
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II  ihc  passage  opens  instantaneously,  the  flowtield  in  the 
passaye  will  have  only  one  discontinuity  —  the  interlace  be- 
tween the  fresh  air  and  exhaust  gas  entering  and  leaving  the 
passaye,  respectively  To  model  this  condition  lor  the  gradual 
passaye  opening,  it  was  assumed  thai  in  the  passage:  /'„  =  I 
aim,  f>„  0.5  kg/in',  U,t  -  I  50  m/s,  and  V„  --  0;  in  the  port: 
/',  I  aim,  ,>,,■={. 4  kg/in',  {Jit-  Is"  m/s,  and  Vlt  0.  The 
velocity  ol  the  passage  opening  was  V  —  200  ni/s. 

Kcsulis  lor  tins  simulation  are  presented  in  fig.  8    It  can  be 
seen  thai,  since  at  first  moment  most  ol  the  inlet  cross  section 


is  blocked  by  the  wall,  a  rarefaction  wave  reflects  from 
inlet  wall.  The  passage  lull  opening  occurs  at  /  =  0. 1  ms.  |- 
this  time  on  the  pressure  deviation  in  the  flowtield  wea 
and  at  /  =  0.22  ms  the  pressure  is  almost  eompletelv  uuil 
and  equal  to  the  static  pressure  in  both  Ihe  undislui  I  v.  I  Ji 
air  ami  exhaust  gas  I  igtiie  Sb  shows  that  the  interlace 
carry  Ihe  imprint  ol  the  gradual  passage  opening  a  long 
alter  the  passage  becomes  fully  open.  There  is  no  dissip; 
mechanism  included  in  the  iuaihematic.il  model  used  in 
study  to  force  the  interlace  to  become  normal  to  the  pas 
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walls.  Therefore,  the  interface  will  "remember"  the  dynamics 
of  the  gradual  passage  opening  in  passages  with  large  length- 
to-width  ratio. 

Conclusions 

It  is  shown  in  the  present  study,  on  the  basis  of  numerical 
modeling,  that  the  dynamics  of  the  passage  opening  sig- 
nificantly affects  the  flow  pattern  in  the  passages  of  wave 
rotor  devices.  For  rectangular  axial  passages,  even  when  the 
velocity  of  the  passage  opening  is  500  m/s,  a  one-dimensional 
flow  pattern  forms  only  in  passages  with  a  length-to-width 
ratio  larger  than  3.  In  the  region  preceding  the  formation  of 
the  one-dimensional  flow  pattern,  the  flow  is  rotational  and, 
in  some  instances,  contains  shock  and  pressures  waves  which 
repetitively  reflect  from  the  upper  and  lower  walls  of  the 
passage.  In  practice  this  would  lead  to  significant  mixing 
between  the  driver  gas  (e.g.,  exhaust  gas)  and  driven  gas  (e.g., 
fresh  air)  and  reduce  the  efficiency  of  the  engine  cycle.  With  a 
reduction  in  the  velocity  of  the  passage  opening,  the  volume 
of  the  mixing  region  increases  as  a  result  of  the  rotational 
flow  which  develops  from  the  slow  opening  of  the  passage. 
With  an  increase  in  shock  wave  strength,  the  volume  of  the 
mixing  region  increases  as  a  result  of  the  rotational  flow 
which  develops  form  multiple  reflections  of  the  shock  and 
pressure  waves  from  the  passage  walls. 

Numerical  modeling  of  the  process  of  gradual  passage 
opening  at  the  inlet  port  of  the  driven  gas  revealed  that  the 
interface  between  the  gases  will  move  all  the  way  through  the 
passage  with  a  frozen  pattern  of  distortion  or  obliqueness  to 
the  passage  walls.  The  interface  obliqueness  increases  when 
the  velocity  of  the  passage  opening  decreases. 


In  all,  it  can  be  concluded  that  taking  into  account 
gradual  passage  opening  is  essential  for  the  wa\e  mac 
design,  not  only  for  proper  timing  and  wave  arrangement 
also  because  o(  the  losses  which  will  occur  due  to  mixing 
wave  reflections.  The  gap  between  projected  and  achii 
efficiencies  of  wave  machlnes:'h,  could  be  partialis  M 
neglect  of  the  effects  of  the  gradual  passage  opening  wind 
studied  herein. 
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'Gradual'  bpehihg'  of'  Skewed'  Passages'  in'  Wave'  Motors 


Introduction 

The  problem  of  gradual  opening  of  rectangular  axial  passage  in  wave 
rotors  was  studied  in  our  article1  published  in  the  January  issue  of  Journal 
of  Propulsion  and  Power.   There  we  have  defined  the  problem  of  gradual 
opening  in  the  wave  rotor  passage  and  the  mathematical  model  which  was  used 
to  simulate  the  opening  process.   In  this  article  we  will  assume  that  the 
reader  is  familiar  with  the  reference1 ,  and  it  could  be  regarded  as  an 
extension  to  the  that  reference. 

If  the  wave  rotor  is  used  to  produce  shaft  power,  it's  passages  should 
be  skewed  in  one  form  or  another.   In  this  paper  we  will  analyze  the  gradual 
opening  of  a  skewed  passage  and  will  examine  the  conclusion  drawn  for  the 
rectangular  axial  passage  for  this  more  general  geometry  of  the  passage. 

Fesults'  and'  biscussibh 

The  main  conclusion  of  the  study  on  gradual  opening  of  rectangular 
passage  is  that  in  order  to  minimize  the  mixing  losses  caused  by  rotational 
flow  in  the  passage,  the  opening  velocity  should  be  very  high.   In  the 
limit,  instantaneous  opening  of  the  wave  rotor  passage  will  lead  to  one 
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dimentional  flow  pattern  in  the  passaoe  will  have  minimal  mixina  losses. 
When  the  passaae  of  the  wave  rotor  is  skewed,  even  instantaneous  openina  of 
the  passage  will'  hbfc  lead  to  development  of  the  one  dimensional  flow  pattern 
with  small  mixina  losses.   That  statement  will  be  demonstrated  in  the 
following  example. 

Let's  model  the  openina  process  for  the  passaae  0.02  m  wide  and  0.24  m 

lona.   It  has  left  and  riaht  hand  inlets  parallel  to  the  v-axis  and  upper 

and  lower  walls  of  the  passaae  form  60  anale  with  positive  direction  of  the 

x-axis.   It  is  assumed  that  initially  air  in  the  passaae  is  at  the  followina 

conditions. 

P   =  1  atm;  P   =1.2  ka/M3;  U   =  0;  V   =0 
o  o        '  o       o 

The  driver  aas  enterina  throuah  the  port  at  the  left  hand  end  was  assumed  to 

have  the  followina  properties: 

P,  =  1.8  atm;  P,  =  1.P1  ka/M3;  n  =  75  M/sec;  v.  =  129.9  M/sec 
d  d  d  d 

The  conditions  for  the  driver  and  driven  aas  are  the  same  as  in  case  of  the 
rectanoular  passaae  which  we  reported  in  the  previous  paraaraph. 

In  Fioures  1a  and  1b  results  for  simulation  of  the  instantaneous 
openina  of  the  skewed  passaae  are  shown  in  form  of  pressure  and  velocity 
contours  at  the  seauences  of  times.   The  flow  pattern  near  the  inlet  in 
Fioure  1  is  hiahly  rotational  which  suaaests  very  hiah  mixina  losses.   That 
is  caused  partially  by  reorientation  of  the  shock  wave.   At  the  first 
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moments  of  time  after  the  openina  of  the  passaae  the  shock  wave  between  the 
driven  and  driver  aases  is  obligue  to  the  lower  and  upper  walls  of  the 
passaae.   At  later  times  this  shock  wave  turns  and  becomes  normal  to  the 
upper  and  lower  walls.   Thus,  incontrast  with  rectanqular  for  skewed 
geometry  passages  there  is  no  obvious  condition  for  minimizing  the  mixing 
losses  caused  by  the  inlet  opening. 

Let's  find  the  conditions  for  opening  of  the  skewed  passaae  which  will 
lead  to  minimal  mixing  between  driver  and  driven  gases.   As  we  have  stated 
before,  rotational  flow  in  the  instantaneously  opened  skewed  passage,  will 
develop  because  of  the  rotation  of  the  shock  wave  from  the  parallel  position 
to  the  inlet  (initially)  to  the  position  where  it  is  normal  to  the  upper  and 
lower  walls  of  the  passage  (when  the  flow  is  developed).   If  we  could  form  a 
shock  wave  which  will  be  all  the  time  normal  to  the  lower  and  upper  passaae 
walls  then  the  rotational  flow  and  mixina  will  be  minimal. 

Analysina  the  formation  of  the  shock  wave  front  at  the  inlet  we 
concluded  that  for  the  shock  wave  formina  at  the  inlet  to  remain  normal  to 
the  lower  wall,  inlet  should  be  openina  at  the  rate: 

v   »  V  .  /sin  a  (1  ) 

op    sh  '     sk 

where  V   -  velocitv  of  the  inlet  openina. 
op 

V   -  shock  wave  velocitv  in  the  media, 
sn 

a    -  anale  of  the  skewed  passaae. 
The  velocitv  is  eaual  to  the  velocitv  with  which  the  shock  wave  surface  will 
slide  alona  the  inlets  wall. 

In  Fias .  2a  and  2b  results  for  the  modelina  of  the  oradual  openina  of 
the  skewed  passaae  inlet  with  the  openina  velocity  in  accordance  to  eauation 
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(2)  are  shown.   As  it  is  obvious  from  coutour  plots  of  pressure  and  velocity 

in  Fiqs.  2a  and  2b,  opening  the  passaoe  with  velocity  V   calculated  from 

op 

equation  (1)  assures  development  of  the  flow  pattern  with  minimal  mixina. 
The  shock  wave  surface  in  that  case  remains  at  all  times  normal  to  the  walls 
of  the  passage. 

&bhclusibhs 

Modeling  of  the  gradual  opening  of  the  skewed  passage  revealed  the  way 
to  minimize  mixing  losses  at  the  passage  inlet.  The  mixing  will  be  minimal 
when  the  velocity  of  the  opening  is  matched  with  shock  wave  velocity  in  the 
passage  divided  by  the  angle  of  skewed  passage.  Oar  simulations  shows  that 
even  when  conditions  for  the  optimal  opening  are  not  satisfied  exactly, 
reduction  of  the  rotational  mixing  could  be  significant  if  the  opening 
velocity  is  +15%  of  the  optimal  value. 

We  can  also  conclude  that  a  very  high  opening  velocity  is  reguired  in 
order  to  obtain  low  mixing  loss.   The  opening  velocity  of  the  passacje  should 
be  always  higher  than  the  velocity  of  the  shock  wave  generated  in  the  wave 
rotor  passage  at  the  high  pressure  inlet.   Fven  substantially  skewed  passaae 
will  reouire  the  opening  speed  °10%  hiaher  than  the  shock  wave  velocity  in 
the  passage  which  is  not  easy  to  achieve  for  the  typical  flow  conditions  in 
the  wave  rotor. 

Another  limitation  is  that  even  when  this  optimal  speed  of  openina  is 
achieved  at  one  port,  it  will  not  be  optimal  at  other  ports;  so, 
optimization  will  not  be  full.   However,  mixing  losses  in  the  skewed 
passages  will  be  smaller  than  in  rectangular  one,  if  the  gradual  opening  of 
the  skewed  passage  begins  from  its  acute  angle. 
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APPENDIX  B 

Numerical  Modeling  of  the  Non-Steady  Thrust  Produced  by  Intermittent  Pressure 
Rise  in  a  Diverging  Channel. 
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ABSTRACT 

The  dynamics  of  the  expansion  of  detonation 
products  in  a  diverging  channel  were  investigated 
numerically.   A  two-dimensional  unsteady  Euler  code 
based  on  the  Godunov  method  was  employed.   The 
influences  of  particular  features  in  the  expansion 
process,  such  as  the  presence  of  reversed  flow  and 
propagation  of  hammer  shocks,  on  the  production  of 
thrust  were  examined.   Sequentially  expansions  of 
detonation  products  were  also  modeled  and  it  was 
concluded  that  in  order  to  maintain  a  high  frequency 
periodic  mode  of  operation  for  propulsion  applications 
the  channel  should  be  refilled  with  ambient  air  after 
each  expansion.  The  influence  of  the  ratio  of  ambient 
air  to  detonation  product  volume  on  the  dynamics  of  the 
thrust  production,  and  on  the  impulse  generated  during 
the  expansion,  are  also  reported. 

NOMENCLATURE 

dAx  -  element  of  channel  internal  surface  area  normal 

to  the  x-axis 
e   -  energy  per  unit  volume 
p   -  pressure 
t   -  time 

u   -  component  of  velocity  in  x  direction 
v   -  component  of  velocity  in  y  direction 
x   -  cartesian  coordinate 
y   -  cartesian  coordinate 
e   -  internal  energy  per  unit  mass 
Y   -  ratio  of  specific  heats 
p   -  density 
subscripts 

0   -  conditions  at  t=0  in  combustion  products 
ex  -  external 
L   -  along  inner  wall  of  the  channel 

INTRODUCTION 

Early  in  the  development  of  propulsion  engines  for 
aircraft,  a  choice  had  to  be  made  between  engine 
concepts  based  on  nonsteady  or  on  steady  gasdynamic 
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processes.   Concepts  based  on  steady  gasdynamic 
processes  were  pursued  largely  because  they  were 
simpler  to  analyze  and  to  appreciate  conceptually.   On 
the  basis  of  the  chosen  steady  concepts,  engines  were 
built,  developed  and  perfected  with  time.   The 
nonsteady  concepts,  which  looked  as  promising 
initially,  were  not  developed  and  remain  today  in  the 
conceptual  stage. 

The  need  for  a  more  thermally  efficient  small 
engines  for  various  military  and  civilian 
applications,  has  led,  at  intervals,  to  a 
reexamination  of  nonsteady  concepts  (1,2,3).   Today, 
such  an  examination  can  be  made  more  thoroughly.   For, 
while  nonsteady  engine  concepts  were  at  one  time 
extremely  difficult  to  analyze,  numerical  modeling  on 
computers  can  now  provide  time-dependent  pictures  of 
internal  flow  processes.   This  can  provide  efficient 
tools  for  preliminary  design  calculations  and 
eventually  analysis  tools  for  design  optimization. 

From  the  late  50' s  until  the  early  70' s,  the 
feasibility  of  an  engine  operating  with  intermittent 
detonative  combustion  was  studied  at  the  University  of 
Michigan  (4,5 ,6) .   Specific  impulses  over  2100  sec. 
were  realized  for  a  single  linear  shock  tube  operating 
intermittently  with  frequencies  up  to  35  detonations 
per  second.   Nicholls  et  al.  showed  that  an  engine 
operating  intermittently  on  detonation  waves  will  have 
some  advantages  over  an  engine  operating  on 
deflagration.   In  general,  it  will  have  very  high 
specific  thrust.  The  engine  will  be  very  simple 
mechanically,  and  not  require  precompression  of  the 
mixture  for  efficient  combustion. 

One  of  the  disadvantages  of  the  engine  is  that  jet 
velocities  developed  after  detonative  combustion  are 
very  high.   Efficient  operation  of  the  engine  for 
propulsion  could  be  realized  only  at  high  supersonic 
vehicle  velocities  (Mach  number  of  about  4).   In  order 
to  reduce  the  velocity  of  the  out-coming  jet,  it  was 
proposed  to  use  spinning  detonation  (6^).   But  the 
spinning  detonation  process  proved  to  be  unsuitable 
for  use  in  an  engine  because  it  was  unstable.   Back  et 
al.  (8)    studied  the  feasibility  of  using  detonative 
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THE  MATHEMATICAL  MODEL  AND  NUMERICAL  SOLUTION 


It  is  assumed  that  the  jet-propulsion  nozzles  of 
an  engine  using  detonative  combustion  can  be 
constructed  from  a  number  of  straight  diverging 
channels  which  are  grouped,  or  clustered  together.   In 
one  cycle  of  the  engine,  a  small  volume  of  each  tube 
is  filled  with  the  combustible  mixture  and  undergoes 
detonation.   Then  the  detonation  products  expand  into 
the  stationary  gas  which  fills  the  rest  of  the  tube  at 
ambient  pressure  and  temperature.  The  tube  is 
diverging  and  so  the  shock  wave  decays  in  strength. 
Finally,  a  weak  shock  wave  leaves  the  tube,  and  a 
subsonic  flow  of  gas  discharges  from  the  exit  of  the 
channel. 

It  is  assumed  that  the  flow  is  periodic  along  the 
lines  which  are  a  continuation  of  the  walls  of  the 
channel.   This  implies  that  we  are  modeling  the  flow 
in  a  centrally  located  channel  of  the  cluster. 

It  is  further  assumed  that  the  flow  is  inviscid. 

The  unsteady  two-dimensional  Euler  equations  can 
be  written  in  conservation  law  form  as: 
3U 
3t 


where 


P 

pu 

pv 

e 


F  = 


3F        3G 
3X        3Y 
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pu 

p  +   pu^ 

puv 

• 

(e   +   p)u 

pv 

puv 

p  + 

pv^ 

(e   + 

p)v 
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Here  p  is  the  density,  u  and  v  are  the  velocity 
components  in  the  X  and  Y  coordinate  directions,  p  is 
the  pressure  and  y  is  the  ratio  of  specific  heats. 
The  energy  per  unit  of  volume,  e,  is  defined  as 

2    2 

f           u  +  v  . 
e  -  p(  e  + ^ )  » 

where  e  =  -, rr—  is  the  internal  energy  per  unit  mass 

(y-Op 

We  look  for  the  solution  of  the  system  of  equations 
represented  by  Eq  (1)  in  the  computational  domain  as 
shown  in  Fig  1  in  time  t,  with  the  following 
conditions  at  the  domain  boundaries: 
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a)  Solid  wall  along  segment  1-2,  1-3  and  2-4 

The  condition  on  the  surface  is  defined  by  solving 
the  Riemann  problem  between  the  point  nearest  to 
the  wall  in  the  domain  of  integration  and  its 
mirror  image  in  the  direction  normal  to  the  wall. 

b)  Periodicity  between  segments  3-5  and  4-6 

By  virtue  of  the  central  location  of  the  channel 
flow  at  each  grid  point  along  3-5  is  the  same  as 
at  corresponding  grid  points  along  4-6  in  Fig  1. 


c)   Outflow  along  segment  5-6 

If  the  outflow  is  subsonic  at  the  downstream 
boundary  we  define  only  the  pressure  pout»  and 
for  all  other  flow  parameters  apply  the 
continuation  condition.   That  means  that  uout, 

vout  an<*  pout  are  set  e<lual  to  the  values  of  u,  v 

and  p  one  point  ahead  of  the  downstream 

boundary. 

If  the  outflow  is  supersonic  the  continuation 

condition  is  applied  to  all  parameters  at  the 

downstream  boundary. 

It  is  assumed  that  initially  at  the  time  t=0,  the 
channel  is  filled  with  detonation  products  from 
segment  1-2  to  l'-2'.   Outside  this  region  of  the 
computational  domain,  the  gas  is  assumed  to  be  at 
ambient  conditions. 

The  Godunov  method  has  been  used  to  obtain  a 
numerical  solution  of  the  Eq  (1)  with  the 
described  boundary  and  initial  conditions.   Details  of 
the  implementation  of  the  method  and  boundary 
conditions  are  given  in  (7) . 

An  orthogonal  grid  is  constructed  which  covers 
the  computational  domain  as  shown  on  Fig  1.   It  should 
be  noted  that  although  the  formulation  of  the  problem 
and  its  solution  is  two-dimensional,  with  the  boundary 
conditions  described  above  the  flow  will  be 
essentially  one  dimensional. 

It  is  assumed  that  at  t*0,  the  detonation  wave 
has  passed  through  the  combustible  mixture  and  the 
products  of  combustion  have  uniform  properties.   Thus 
the  time  of  initiation  and  propagation  of  the 
detonation  wave  through  the  combustible  mixture  is 
neglected  in  comparison  with  the  time  for  propagation 
through  the  channel  and  subsequent  expansion  of  the 
gases. 
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RESULTS  AND  DISCUSSION 

The  history  of  the  flow  development  in  the 
diverging  channel  is  shown  in  Fig  2.   Here  the 
velocity  and  pressure  distributions  at  the  center-line 
of  the  channel  are  shown  at  specific  moments  of  time, 
where  time  t=0  is  when  the  expansion  begins.   At  time 
t=0,  approximately  8.5%  of  the  channel  volume  from  the 
left-hand  side,  is  filled  with  detonation  products 
with  initial  pressure  density  and  velocity  having 
values  of 

p0  a  40  atm,  p0  ■  5.2kg/m3,  v0  ■  0  and  u0  =  0 
respectively.  At  t=0,  the  rest  of  the  channel  is 
filled  with  ambient  air  with 

p  =  1  atm,  p  =  1.3  kS/nr^ ,  v  ■  0  and  u  ■  0  . 
The  length  of  the  channel  is  1  meter.  These 
conditions  will  be  referred  to  as  Case  1. 

Initially,  a  shock  wave  going  from  left  to  right 
and  a  rarefraction  wave  going  from  right  to  left  can 
be  seen  in  Fig  2.   The  shock  wave  progressively  decays 
because  the  channel  is  diverging.   The  rarefraction 
wave  propagates  rapidly  towards  the  wall  at  the  left 
end  of  the  channel,  because  the  temperatures  in  this 
region  are  very  high  (approximately  2800  K)  and  the 
channel  is  converging  towards  the  left  end.   At  a  time 
of  about  0.15  msec  the  rarefraction  reaches  the  wall 
at  the  left  end  of  the  channel.  From  this  moment, 
pressure  at  the  left  end  of  the  channel  rapidly 
decreases  from  about  40  atm  at  t  =  0.15  msec  to  about 
0.3  atm  at  t  =  0.82  msec.   Because  of  this  rapid 
depressurization  of  the  region  adjacent  to  the  left 
end  wall,  the  pressure  on  the  right  hand  end  of  the 
channel  becomes  higher  than  that  at  the  left,  which 
generates  a  shock  wave  going  from  right  to  left  at 

time  fO.fi  msec.   This  back-going  shock  wave  rapidly 
slows  down  and  then  reverses  the  direction  of  the  flow 
because  the  channel  is  converging  towards  its  left  end 
wall. 

The  negative  flow  reflects  from  the  left  end  wall 
at  time  t«1.4  msec,  forming  a  hammer  shock  wave  which 
moves  from  left  to  right  and  again  reverses  the  flow 
direction.   The  formation  of  the  hammer  shock  wave  can 
be  seen  clearly  on  the  velocity  and  pressure  graphs 
for  time  t  =  1.6  msec.   Because  the  shock  travels  in  a 
diverging  channel,  it  weakens  rapidly. 

The  main  shock  wave  produced  by  the  detonation 
products  leaves  the  channel  at  the  time  «0.94  msec, 
and  because  the  channel  is  diverging  in  the  direction 
of  propagation  the  pressure  behind  the  shock  wave 
drops  from  about  10  atm  at  time  t«»0.1  msec  to  about  5 
atm,  at  t  =  0.94  msec  when  it  leaves  the  channel. 
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The  mass  velocity  of  the  gas  drops  from  about 
1100  m/sec  at  the  beginning  of  the  expansion  to  about 
500  m/sec  when  the  gas  begins  to  leave  the  detonation 
tube.   The  velocity  of  the  gas  leaving  the  tube  then 
decreases  to  -20  m/sec  and  increases  again  to  about 
50  m/sec  when  the  shock  wave  reflected  from  the  left 
end  wall  reaches  the  exit  of  the  channel.   This 
concludes  the  significant  events  which  occur  in  the 
diverging  channel  as  a  result  of  the  expansion  of  the 
detonation  products.   At  later  times,  rarefraction  and 
pressure  waves  continue  to  form  and  travel  in  the 
channel,  but  they  are  much  weaker  and  their  influence 
on  the  thrust  produced  is  very  small. 

In  Fig  3  thrust  as  function  of  time  is  plotted 
for  the  case  described  in  the  graphs  given  in  Fig  2. 
Thrust  was  calculated  for  the  case  of  a  stationary 
tube  with  external  pressure  equal  to  the  ambient 
pressure  using 

It  can  be  seen  in  Fig  3  that  most  of  the  thrust 
is  produced  in  the  2  msec  following  the  beginning  of 
the  expansion  of  the  detonation  products.   Integration 
of  thrust  in  time  shows  that  90%  of  the  impulse  is 
produced  in  the  first  2  msec  and  10%  in  the  following 
7  msec. 

The  peculiarities  of  the  thrust  curve  can  be 
understood  in  relation  to  the  wave  pattern  in  the 
tube.   The  first  change  in  slope  at  time  t»0.8  msec 
corresponds  to  the  tine  when  the  front  of  the  main 
shock  wave  reaches  the  end  of  the  tube.   Because  the 
pressure  behind  the  shock  front  drops  rapidly,  the 
thrust  decay  increases  when  the  main  shock  front 
leaves  the  channel.   At  time  t«1.5  msec  the  thrust 
increases  as  a  result  of  the  hammer  shock  wave 
reflecting  from  the  left  wall. 

When  the  main  process  of  expansion  ends,  at 
time  t«6  msec,  the  pressure  in  the  tube  is  about  1.1 
atm  and  the  average  density  is  about  0.43  kg/m^.   This 
corresponds  to  an  average  temperature  in  the  channel 
of  about  935  K.   Since  the  pressure  in  the  tube  is 
close  to  ambient,  the  heat  exchange  will  be  by  natural 
convection,  and  it  will  take  a  relatively  long  time 
for  the  gas  remaining  in  the  tube  after  the  expansion 
to  cool. 

In  order  to  examine  what  will  be  the  wave  pattern 
and  thrust  when  a  second  detonation  wave  expands  in 
the  channel  immediately  following  the  expansion  of  the 
first  detonation  wave  (6  msec  after  the  beginning  of 
the  process  of  expansion  of  the  first  detonation 
wave),  conditions  were  simulated  numerically  in  the 
following  way.   Results  of  the  calculation  for  the 
first  detonation  wave  expansion  were  used  as  the 
Initial  distributions  of  flow  parameters.   Then,  in 
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the  8.5%  of  the  channel  volume  adjacent  the  left  end 
wall  the  pressure,  density  and  velocity  were  given  the 
values 

pQ  =  40  atm,  p0  =  5.2  kg/m3  and  v0  ■  u0  »  0 
respectively  at  time  t=0.  These  conditions  will  be 
referred  to  as  Case  2. 

In  Fig  4,  the  thrust  vs  time  is  shown  for  the 
second  detonation  wave  expansion.   Since  the 
detonation  products  expand  in  a  low  density  hot  gas, 
the  expansion  goes  much  faster  than  in  the  previous 
case.  At  time  t»0.8  msec  from  the  beginning  of  the 
second  detonation  wave  expansion,  the  thrust  drops  to 
zero.   This  can  be  compared  with  a  thrust  of 
approximately  6»10^N/m  at  the  corresponding  time  shown 
in  Fig  3,  for  expansion  into  ambient  air.   Then, 
because  the  large  volume  of  the  channel  has  static 
pressure  lower  than  ambient,  thrust  becomes  negative 
and,  as  in  the  previous  case,  a  back-going  shock  wave 
develops.   In  Fig  5,  the  development  in  time  of  this 
recursive  shock  can  be  followed.   It  is  observed  that 
the  negative  flow  velocities  which  develop  are  twice 
as  large  as  those  for  the  preceding  expansion  into 
ambient  air.  And  at  time  t»2  msec  the  flow  throughout 
the  tube  is  reversed.   Reflection  of  the 
left-traveling  shock  wave  from  the  left  end  wall  at 
time  t»2  msec  the  flow  throughout  the  tube  is 
reversed.   Reflection  of  the  left-traveling  shock 
wave  from  the  left  end  wall  at  time  t«2  msec  produces 
a  very  strong  hammer  shock  which  reverses  the  thrust 
and  the  flow  direction  at  time  t«2.1  msec.   Because  of 
the  large  negative  thrust,  the  total  impulse  which  is 
produced  in  Case  2  Is  about  15%  smaller  than  the 
impulse  of  the  detonation  products  expanding  in 
ambient  air. 

Another  notable  effect  of  expansion  into  the 
products  of  the  previous  cycle,  is  the  substantial 
rise  of  temperature  which  occurs  in  the  channel  after 
the  second  expansion.   When  the  process  of  the  second 
expansion  ends  the  average  temperature  in  the  channel 
is  approximately  1550  K.  This  is  615  K  higher  than 
the  final  temperature  after  expansion  into  ambient 
air. 

In  order  to  examine  how  thrust  and  total  impulse 
are  influenced  by  the  volume  of  air  contained  at 
ambient  conditions  in  the  channel  before  the 
expansion,  two  additional  test  cases  were  modeled.   In 
Case  3,  the  volume  of  ambient  air  in  the  channel  to 
the  right  side  of  the  interface  with  the  detonation 
products  is  taken  to  be  approximately  2.5  times  larger 
than  the  volume  of  the  detonation  products.   This 
requires  a  channel  about  50%  shorter  than  the  original 
channel,  for  which  the  ratio  of  the  volumes  of  ambient 
air  to  detonation  products  was  about  11.   In  Case  4, 


there  was  no  air  in  the  channel  when  expansion  began, 
which  implies  that  the  channel  in  this  case  has 
only  about  25%  of  the  original  length.  Therefore, 
between  Cases  3  and  4  and  Case  1  only  the  length  of 
the  channel  was  changed.   The  geometry  of  the  channel, 
initial  parameters  of  the  detonation  products  were  the 
same  in  Cases  1,  3  and  4. 

In  Fig  6  the  thrust  is  shown  as  a  function  of 
time  for  Cases  1,  3  and  4.   It  can  be  seen  that  Case  1 
(solid  line)  where  the  channel  is  fully  extended, 
gives  consistently  higher  thrust  and  that  the  thrust 
decreases  when  the  volume  of  ambient  air  in  the 
channel  (or  channel  length)  decreases.   Case  3  (chain- 
dotted  line)  and  Case  4  (dashed  line)  have  regions 
where  thrust  is  negative.   Comparison  of  the  total 
impulse  generated  shows  that  the  extension  of  the 
channel  is  very  beneficial  for  producing  impulse.   The 

impulse  produced  in  Case  4,  where  the  tube  is  filled 
with  detonation  products  is  0.162«l(pN/M»sec.   When 
the  channel  is  extended  so  that  only  28.6%  of  the 
channel  volume  is  initially  filled  with  detonation 
products  (Case  3),  the  impulse  increases  28%.   An 
additional  extension  of  the  channel  so  that  only  about 
8.5%  of  the  tube  volume  is  initially  filled  with 
detonation  products  (Case  1),  leads  to  total  impulse 
of  0.2784»103  N/M'sec,  which  is  a  72%  increase 
over  Case  4. 

CONCLUSIONS 

Numerical  modeling  of  the  process  of  expansion  of 
detonation  products  in  a  diverging  channel  revealed  a 
flow  pattern  rapidly  changing  in  time  with  multiple 
shock  and  rarefraction  waves.   One  of  the  most 
interesting  features  was  the  occurrence  of  a 
backwards-traveling  shock  wave,  which  developed  as  a 
result  of  over-expansion  in  the  channel.  The  shock 
wave  first  reversed  the  flow  direction  (and  in  some 
cases  developed  a  negative  thrust)  and  then, 
reflecting  from  the  left  end  wall  as  a  hammer  shock, 
reversed  the  flow  direction  again,  increasing  the 
thrust. 

Modeling  of  the  second  detonation  wave  expanding 
into  the  products  of  the  previous  expansion,  showed 
that  this  arrangement  leads  to:   a)  significant 
negative  thrust;  b)  increases  of  the  gas  temperature 
after  expansion;  c)  losses  of  15%  of  the  total 
impulse.   This  leads  to  the  conclusion  that  the 
channel  should  be  filled  with  ambient  air  after  each 
detonation  and  expansion.   Since  the  process  of 
expansion  takes  only  about  0.006  sec,  the 
time-averaged  thrust  which  can  be  generated  by  a 
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single  channel  will  be  limited  by  the  rate  at  which 
the  channel  can  be  re-charged  with  air  and  mixture, 
and  refired.   For  very  low  initial  velocities  in  the 
channel  (consistent  with  the  initial  conditions 
specified  in  the  modelling),  the  frequency  of  firing 
is  limited  to  about  20  detonations  per  second.   This 
corresponds  to  0.044  sec  for  recharging  at  gas 
velocities  of  about  23  m/sec.   Investigation  of  the 
influence  of  the  channel  length  on  the  impulse  and 
thrust  production  showed  that  increases  in  the  volume 
of  ambient  air  in  the  channel  ahead  of  the  detonation 
products,  which  was  obtained  with  an  increase  in 
channel  length,  led  to  significant  increases  in  the 
thrust  and  impulse  generated  during  the  expansion. 
The  total  impulse  generated  by  an  expansion  in  the 
channel  in  which  8.5%  of  the  total  volume  was  filled 
initially  by  the  detonation  products,  was  found  to  be 
72%  greater  than  in  the  case  of  a  short  channel  filled 
only  with  the  same  volume  of  detonation  products  at 
the  same  initial  conditions. 

The  diverging  channel  geometry  provided  the 
simplest  geometry  with  area  change  to  model  using  the 
unsteady  2D  Godunov-Euler  code.   Variations  in  the 
geometry  can  now  be  made  in  order  to  examine  the 
unsteady  thrust  performance. 
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Figure  1.   The  Computational  Domain. 
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Figure  2.   The  Evolution  in  Time  of  the  Detonation 
Products  Expansion. 
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Figure  3.   The  Evolution  in  Time  of  the  Thrust 

Developed  by  the  Expansion  of  the  Detonation 
Products. 
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Figure  4.   The  Thrust  vs.  Time  Curve  for  the  Second 
Sequential  Expansion  of  the  Detonation 
Products. 
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Figure  5.   Formation  and  Evolution  of  the  Negative 

Flow  Velocities  for  the  Second  Sequential 
Expansion  of  the  Detonation  Products. 
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Figure  6.   The  Thrust  vs.  Time  Curves  for  Expansion 
in  Channels  of  Different  Lengths. 
Solid  Line  -  Longest  Channel 
Dashed  Line  -  Shortest  Channel 
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INTRODUCTION 

Wave  Rotors  are  devices  which  use  unsteady 
pressure  waves  to  effect  direct  energy  exchange  between 
two  gases.   They  offer  the  potential  for  significant 
Improvements  in  air-breathing  engine  propulsion  cycles 
through  their  capability  of  withstanding  higher 
temperatures  and  pressures  than-  present-day  turbines. 
Other  diverse  applications  are  recounted  in  Ref  (1;. 

In  its  most  basic  form,  a  wave  rotor  is  a  drum 
with  axial  or  helical  (staggered)  passages  arranged 
around  the  periphery.   This  single  drum-like  rotor 
replaces  separate  compressor  and  turbine  components  for 
gas  turbine  applications.   The  compression  and 
expansion  of  the  two  fluids  occurs  in  the  passages  as  a 
result  of  shock  tube  like  processes.   In  a  typical 
configuration,  combustion  products  at  high  temperature 
and  pressure  give  up  energy  to  the  'driven'  fluid 
(usually  air)  through  the  action  of  time-unsteady 
compress ion /shock  waves.   The  combustion  products,  (hot 
'driver'  gases),  are  in  turn  expanded  and  exhausted 
from  the  rotor,  the  available  work  of  expansion  being 
utilized  to  Induce  a  fresh  charge  of  air  onto  the 
rotor.   Careful  sequencing  of  the  passage  ends  past 
stationary  inlet  and  outlet  ports  in  valve  plates  on 
either  side  of  the  rotor  creates  a  cyclic  Internal  wave 
pattern  in  the  wave  rotor  component.   The  high 
temperature  capability  of  the  device  Is  a  direct 
consequence  of  the  repetitive  processing  of  both  cold 
and  hot  fluid  alternately  in  the  same  rotor. 

Estimation  of  the  aero- thermodynamic  performance  of 
these  devices  hinges  on  calculations  of  the  unsteady 
energy  exchange  in  the  wave  rotor  component.   Numerical 
simulation  of  the  unsteady  wave  processes  can  generally 
be  carried  out  on  a  one-dimensional  basis.   However, 
the  complex  pattern  of  flow  discontinuities  and  wave 
Interactions  known  to  exist  in  wave  rotors  call  for 
numerical  methods  which  solve  the  nonlinear,  hyperbolic 
system  of  governing  equations  without  relying  on  either 
artificial  viscosity  or  special  treatment  of 
discontinuities.   Described  In  this  paper  are  two  such 
techniques,  the  Codunov  method  and  a  Random  Choice 
method  (Gllmm's  method  or  Plecewlse  Sampling  method). 


NUMERICAL  FORMULATION 

An  essential  component  of  both  the  Random  Choice 
method  and  the  Codunov  method  is  the  solution  of  a 
sequence  of  local  Rlemann  problems.   A  Rlemann  problem 
is  defined  by  setting  up  an  initial  value  problem  for 
the  equations  of  motion  in  Eulerian  form.   The 
unsteady,  two-dimensional  Euler  equations  can  be 
written  in  vectorized,  conservation  law  form  as: 

Ut  +  IP(U)]X  +  [G(U)]y  -  0,  where  U  -  [ p  pu  pv  e]T, 

P(U)  -  [pu  (p  +  pu2)  puv  (e  +  p)u]T,  and 

G(U)  -•  [pv  puv  (p  +  pv2)  (e  +  p)v]T  . 

Here,  p  is  the  density,  u  and  v  are  the  velocity  com- 
ponents in  the  x  and  y  coordinate  directions,  and  p  is 
the  pressure.   Subscripts  indicate  partial  differenti- 
ation with  respect  to  that  Independent  variable.   The 
energy  per  unit  volume,  e,  is  defined  by: 

e  -  p(  e  +  (u2  +  v2)/2),  where  e  -  p/(y-l)p 

is  the  Internal  energy  per  unit  mass. 

The  Rlemann  problem  is  intrinsically 
one-dimensional  in  nature  and  is  defined  accordingly. 
For  the  one  space  variable  case,  thus,  U  -  U(X,t),  and 
the  initial  value  problem  is  set  up  by  specifying 
initial  conditions  which  consist  of  Intervals  over 
which  data  are  constant,  separated  by  jump 
discontinuities.   If  the  time  step  is  chosen  to  be 
sufficiently  small,  the  waves  propagating  with  finite 
speeds  from  adjacent  discontinuities  remain  within 
their  respective  spatial  cells  and  do  not  Intersect. 
This  sequence  of  solutions  from  adjacent  Rlemann 
problems  Is  then  pieced  together  to  obtain  the  whole 
solution  at  each  succeeding  time  step. 

The  Random  Choice  method  (RCM)  for  gasdynamics  is 
the  outcome  of  a  constructive  existence  proof  of 
solutions  to  systems  of  nonlinear  hyperbolic 
conservation  laws  presented  by  Glinm  Ref  (2),  and  Its 
development  Into  an  effective  numerical  tool  by  Chorln 
Refs  (3,4).   The  solution  Is  advanced  In  time  by  a 
method  that  Includes  a  solution  of  Rlemann  problems  as 
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described  above  and  a  sampling  procedure,  which 
chooses  values  from  a  representative  point  of  the 
exact  solution  to  the  local  Rlemann  problem.   The 
sampling  procedure  is  due  to  van  der  Corput  Ref  (5) 
and  generates  equldlstributed  sequences  Ref  (6).   The 
technique  eliminates  numerical  diffusion  and  flow 
discontinuities  are  computed  with  infinite  resolution, 
i.e.,  there  is  no  smearing. 

The  Godunov  method  is  a  finite  difference  method 
which  computes  the  Rlemann  problem  solutions  as  a 
first  step,  and  using  the  fluxes  thus  obtained, 
advances  in  time  by  solving  the  first-order  finite 
difference  approximation  of  the  Euler  equations. 

Coding  details  for  the  RCM,  and  the  Godunov 
method  are  given  In  Refs  (4,7)  and  Refs  (7,8) 
respectively.   In  the  following  section,  examples  of 
the  applications  of  these  techniques  to  wave  rotor 
devices  are  given. 

EXAMPLE  APPLICATIONS 

Fig  1  shows  a  wave  diagram  (a  one-dimensional 
space-time  plot)  for  a  wave  rotor  device  configured  to 
operate  as  a  turbine  alone.   Air  at  a  pressure  higher 
than  ambient  enters  the  rotor  through  the  inlet  port, 
generating  a  shock  wave  with  a  slower  moving  contact 
discontinuity  behind  it  at  point  'a'.   The  shock 
reflects  off  a  wall  boundary  at  point  'b',  crosses  the 
incoming  interface  at  point  ' c'  (bringing  it  to  a 
halt),  and  reaches  the  inlet  side  again  at  point  'd', 
whereupon  the  inlet  port  is  closed.   The  gases  in  the 
rotor  passages  are  now  in  an  essentially  quiescent, 
high  pressure  state.  At  point  'e',  the  outlet  port  is 
opened  to  a  lower  pressure,  a  rarefaction  fan  being 
generated  in  the  process.   The  expansion  waves  travel 
to  the  left,  reflect  off  the  solid  boundary  and 
propagate  back  to  the  outlet  side.  The  port  is  dosed 
when  conditions  in  the  passages  are  essentially  the 
same  as  at  the  beginning  of  the  cycle. 
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Pig.  1   Wave  Diagram  Computed  by  1-D  Random  Choice 
Method. 

S-Shock;  RS-Reflected  Shock;  R-Raref action 
Pan;  RR-Reflected  Rarefaction;  I-Interface 


The  entire  wave  diagram  was  generated  by  the 
one-dimensional  sampling  method  with  properly 
implemented  boundary  conditions.   Pig  2  shows  a 
sequence  of  the  propagation  of  the  shock  wave  and 
contact  discontinuity  in  time.  The  perfect  resolutic 
of  the  discontinuities  is  noteworthy,  even  when  a  vei 
small  density  change  occurs  across  the  Interface,  as 
In  this  case. 
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Pig.  2   Sequence  Showing  Shock  and  Interface  Movement 
S  -  Shock;  I  -  Interface;  RS  -  Reflected 
Shock;  N  -  Times tep 

The  transient  process  of  the  rotor  passage  end 
opening  or  closing  creates  losses  which  affect  the 
performance  of  the  device.   The  gradual  (as  opposed  I 
Instantaneous)  opening/closing  of  the  passages  is 
essentially  a  two-dimensional  flow  phenomenon  and  is 
modelled  using  the  2-D  Godunov  code.   Pig  3  shows  a 
sequence  of  pressure  contours  for  the  gradual  opening 
case  which  clearly  shows  the  significant  effect  the 
dynamics  of  the  passage  opening  has  on  the  flow 
pattern  In  wave  rotor  devices.  The  actual  shock 
formation  for  the  modelled  case  is  seen  to  occur 
approximately  halfway  into  the  passage  as  opposed  to 
the  Instantaneous  formation  generally  assumed  for  th< 
Ideal  case. 

DISCUSSION 

The  two  techniques  described  here  are  powerful 
tools  for  the  design  of  wave  rotor  devices.   The 
example  application  of  the  one-dimensional  Random 
Choice  method  deals  with  a  fairly  elementary  wave 
diagram,  but  serves  to  illustrate  Its  potential  and 
suitability  for  designing  complex  gas  turbine  engine 
cycle  type  wave  diagrams.   The  Godunov  method  is  mor« 
amenable  to  extension  to  multidimensional  analysis 
required  for  loss  calculations,  while  fulfilling  the 
requirements  outlined  in  the  introduction. 
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Fig.  3  Contour  Plots  Generated  by  2-D  Godunov  Method 
for  Gradual  Opening  of  Wave  Rotor  Passage  End. 
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APPENDIX  D 

Local  Cell  Orientation  Method  -  Illustration  of  the  Method  for  the  Solution 
with  an  Oblique  Shock. 
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Local    Cell    Orientation    Method. 
Illustration   of  the  method   for  the   solution   with   oblique   shock 


INTRODUCTION 


The  Finite  Volume  Schemes  for  numerical  solution  of  the  Euler  equa- 
tions  can    be   split   into  two  main    steps: 

a)Calculation   of   the  fluxes   on   the  cell   edges; 

b)Updating   the   values   at  the  center  of   the  cell    ,  using   the 
calculated   fluxes,    to   satisfy   the   Euler  conservation    laws. 

Usually  the  accuracy  of  the  flux  calculation  determine  the  accuracy 
of  the  numerical  modeling  Ref.l.  In  the  region  of  the  shock  wave  the 
accuracy  of  the  flux  calculation  is  significantly  reduced  and  in  some 
methods  special  care  should  be  taken  in  order  to  reduce  the  pre-shock 
and  after-shock  oscillations  .When  the  upwind  methods  are  used  the 
fluxes  are  usually  calculated  using  the  one  dimensional  wave  propagation 
problem  (i.e.  Riemann  problem).  In  this  case  the  assumption  of  the 
one-dimensionality  lead's  to  large  computational  errors  and  significant 
shock    smearing    in    the    region   of  the  oblique   shock. 

In  the  current  study  it  will  be  demonstrated  that  the  accuracy  of 
the    numerical    modeling    could    be   significantly    improved    by    proper    Local 
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PAGE  2 
Cell  Orientation  (SELCO  Method).  We  applyed  the  SELCO  Method  only 
with  the  Godunov  scheme  for  the  solution  of  the  Euler  equations,  but 
the  idea  of  the  method  is  general  enough  to  be  used  with  other  finite 
volume  methods. 


THE  SELCO  METHOD 


For  the  Godunov  method  the  numerical  fluxes  are  calculated  using 
the  solution  of  the  Riemann  problem.  First  the  solution  is  assumed  to  be 
piece-wise  constant  in  every  cell.  Then,  to  calculate  the  flux  for  every 
cell  edge  the  Riemann  problem  is  solved  between  the  cells  adjacent  to 
the  edge.  Solution  of  the  Riemann  problem  is  obtained  in  the  direction 
normal  to  the  cell's  edge.  The  equations  are  then  integrated  in  the 
each  cell  for  one  time  step  using  the  finite  difference  approximation  of 
the  Euler  equations  and  the  fluxes  at  the  cell's  edges.  A  more  complete 
description   of  the  Godunov   method   can   be  found    in    Ref.1. 

In  the  two  dimensional  Godunov  method  it  is  assumed  that  the  fluxes 
could  be  determined  using  the  solution  of  the  one  dimensional  Riemann 
problem  if  the  CFL  number  used  is  smaller  than  0.5.  In  the  following 
example  it  will  be  demonstrated  that  that  assumption  lead's  to  consider- 
able  errors   for  the   solutions   with   oblique   shock   waves. 

As  a  test  case  the  supersonic  flow  over  the  22.3°  wedge  was  simu- 
lated. The  grid  and  the  computational  domain  is  shown  on  the  Figure  1. 
The  system  of  the  two  dimensional  Euler  equations  as  in  Ref.1,  was 
solved   by    time  marshing    in    the   computational    domain    shown    in    Fig.l    for 
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time  t*oo  with    the  following   conditions    at  the   domain    boundaries: 

a)  Inflow  of  the  air  at  M=2.2   along    segment   1-2; 

b)  Outflow   along   segment  3-4; 

(the  flow   along   this   segment  will    be   supersonic   and 
the  continuation   condition    is   applied   to   all    parameters 
at  the  downstream   boundary) 

c)  Solid   wall   along   segments    1-3   and   2-4. 

In  the  Figure  2  (a,b,c,d)  results  of  the  flow  simulation  using  the 
Godunov  method  in  the  computational  domain  shown  on  Fig.1  (squares) 
are  compared  with  the  analytical  solution  for  the  supersonic  flow  over 
the  wedge  (solid  lines)  Ref.2.  Comparison  is  done  for  the  data  on  the 
lower  surface  of  the  channel  shown  on  the  Fig.l.  It  can  be  concluded 
from  this  comparison  that  only  pressure  coefficient  is  calculated  accu- 
racy by  the  Godunov  method.  Especially  notable  the  error  in  entropy. 
Exact  entropy  change  on  the  shock  wave  is  2  times  smaller  than  that 
predicted  by  the  Godunov  method.  The  isomach  lines  are  shown  for  the 
same  simulated  case,  in  Figure  2  (e) .  The  shock  is  smeared  over  the 
large  area  of  the  mesh  (up  to  5-6  points),  but  the  shock's  angle  is 
simulated   correctly. 

In  order  to  study  the  source  of  errors  for  this  type  of  problems  the 
supersonic  flow  over  the  wedges  with  different  angles  was  simulated.  It 
was  found  that  the  pressure  was  predicted  accuratly  for  the  all  simu- 
lated cases.  At  the  same  time  the  error  in  entropy  increases  when  the 
shock  wave  obliquity  increases.  Dependence  of  the  accuracy  of  the 
shock   simulation    on    the   shock  obliquity    is    not   surprising    since   the   basic 
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assumption  about  independance  of  the  flux  calculation  for  the  intersect- 
ing edges  of  the  cell,  for  t  <  0.5*CFL,  became  incorrect  in  the  vicinity 
of  the  shock.  Spliting  the  flux  calculation  for  each  cell  on  four  one 
dimensional  Riemann  problems  for  each  of  the  cell's  edges  would  not 
lead  to  the  large  errors  in  the  regions  of  monotonic  flow.  It  will  also 
give  accurate  approximation  in  the  region  of  the  oblique  shock  if  the 
shock  wave  is  parallel  to  one  of  the  cell  edges.  The  last  statement  will 
be   illustrated   by   the  following    numerical   example. 

Let's  solve  the  problem  of  the  supersonic  flow  with  M  =2.2  over  the 
22.3°  wedge,  as  in  previous  example,  on  the  grid  shown  in  Figure  3. 
The  grid  in  Fig.1  is  skewed  uniformly  so  that  vertical  lines  of  the  grid 
are  forming  a  51.7°  angle  with  the  positive  direction  of  the  x-axis.  An 
angle  of  51.7°  corresponds  to  the  shock  wave  angle  calculated  analyt- 
ically Ref.2.  Because  the  grid  is  skewed  the  oblique  shock  waves  will 
be  parallel  to  the  vertical  grid  lines  of  the  mesh  shown  in  Fig. 3. 
Results  of  this  test  case,  in  the  form  of  the  distribution  of  pressure 
coefficient,  entropy,  density  and  velocity,  are  compared  to  the  analyt- 
ically calculated  values  (solid  lines)  in  Figure  4  (a,b,c,d).  In  Fig. 4 
(e),  the  isomach  lines  for  this  simulation  are  shown.  It  can  be  con- 
cluded from  the  results  presented  in  Fig. 4,  that  the  flow  is  modeled  is 
this  case  with  very  high  accuracy.  Shock  wave  is  resolved  on  one  grid 
cell  and  entropy  jump  is  calculated  exactly.  All  that  is  result  of  the 
proper  cell   orientation   towards   the   shock   wave. 

Skewing  the  entire  grid  can  help  to  accuratly  resolve  only  one 
shock  wave,  and  only  on  the  condition  that  the  shock  wave  angle  is 
known    prior    to    the    solution,     so   it    could    not    be    applied    effectivly.     But 
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if  skewing  of  the  cell  can  be  done  locally  in  the  region  of  the  oblique 
shock   it  can    improve  the   accuracy   of  the   shock   wave   simulation. 

The    proposed    SELCO    method   does    the  cell    reorientation    locally.    The 
method   consist  from   the  following   steps: 

1 )  Integration   of  the   Euler  equations   by   the-Godunov   method; 

2)Define  the  approximate   shock   location    and   the   shock   angle   using 
the  expressions   for  oblique   shocks    Ref.2. 
(Based   on   the  fact  that   pressure   is   calculated   accuratly 
by   the  Godunov   method); 

3)  Rotate  the   cell    edges   that   are   directed   along   the   shock   about 
their  middle   points   on    the   shock   wave   angle.    So,    after    rotation, 
two  of   the   cell    edges   will    be   parallel    locally   to   the   shock   wave 
surface    (see    Figure  5); 

4)Calculate  the  fluxes   on   the   new   edges   of   the  cell; 

5)lntegrate   the   Euler  equations    in   the   new   cell. 

All  these  additional  steps  do  not  add  much  of  computational  work, 
because  the  cell  reorientation  should  be  done  only  in  the  vicinity  of  the 
shock  wave.  If  the  shock  could  be  resolved  on  one  grid  point,  only  one 
cell  would  be  transformed.  Our  experience  has  shown  that  there  is  no 
need  to  extend  the  SELCO  procedure  on  the  cells  ahead  and  behind  the 
shock   wave   surface. 

In  Figure  5  (a,b,c,d)  results  are  shown  for  the  same  flow  condition 
as    in    previous   cases:    supersonic   flow   with    M   =2.2  over   the   wedge  of 
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22.3°.  Calculations  were  done  on  the  grid  shown  in  Figure  1  and  in 
this  case  the  SELCO  method  was  applyed.  It  can  be  concluded  from  the 
results  presented  in  Fig.  5  that  the  accuracy  of  the  shock  wave  model- 
ing using  the  SELCO  method  approach  that  demonstrated  in  Figure  4  for 
the  completely  skewed  grid,  and  is  superior  to  the  accuracy  of  the 
standard  Godunov  method.  The  isomach  lines  for  the  case  where  the 
SELCO  method  was  applied  are  presented  in  Figure  5  (e)  .  The  shock 
tickness  in  this  graph  is  minimal  and  much  more  improved  compared  to 
the   simulation    by   the  original   Godunov   method    shown    in    Figure  2    (e)  . 


CONCLUSIONS 


It  has  been  demonstrated  that  inaccurate  modeling  of  the  oblique 
shock  waves  produced  by  the  Godunov  method  is  the  result  of  the 
obliqueness  of  the  shock  wave  with  respect  to  the  edges  of  the  cells  of 
the  computational  grid  covering  domain  of  integration.  It  is  also  shown 
that  only  when  the  shock  surface  is  parallel  to  the  two  opposite  edges 
of  the  cell   the  oblique   shock   can    be  accuratly   calculated. 

A    new   method   of  the   Local    Cell    Orientation    (SELCO  Method) 
is    proposed     in    order    to    allow    local     reorientation    of    the    cells     in     the 
vicinity    of    the     shock    waves.     The    efficiency    of    the     SELCO    method     is 
demonstrated     for     the     simulation     of     the    oblique     shock     waves     in     the 
supersonic   flow    using    Euler  equations. 

Although   the   new   SELCO   method   was   demonstrated   with   the   Godunov 
scheme    it    will    be    effective    in    applications    to    the    other    upwind    methods 
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that   use  the   finite   volume  formulation. 
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FIGURE  CAPTIONS. 

1.  The   Computational    Domain. 

2.  Computation    by   the   Godunov   method  .  M-Oo  =2.  2. 

Wedge  Angle   is   22.3°.    Continuous    lines 
show  the  exact   solution   for   -    a,b,c   and   d. 

a)  Surface   Pressure   Coefficient. 

b)  Surface   Entropy. 

c)  Surface   Density. 

d)  Surface  Mach    Number. 

e)  Isomach    Lines. 

3.  The   Computational    Domain   with   the   Skewed   Grid. 

4.  Computation   by   the  Godunov   Method   on   the   Skewed   Grid   Shown 

in    Figure  3.    M_oo=2.2.    Wedge   Angle    is   22.3°. 

Continuous    lines    show  the  exact   solution   for   -   a,b,c   and   d. 

a)  Surface   Pressure   Coefficient. 

b)  Surface   Entropy. 

c)  Surface   Density. 

d)  Surface  Mach    Number. 

e)  Isomach    Lines. 

5.  The    Local    Cell    Reorientation    by   the   SELCO   Method. 

6.  Computation    by   the   SELCO   Method  .M_^=2.  2. 

Wedge   Angle   is   22.3°.    Continuous    lines 
show   the   exact   solution   for   -    a,b,c   and   d. 

a)  Surface    Pressure   Coefficient. 

b)  Surface   Entropy. 

c)  Surface   Density. 

d)  Surface  Mach    Number. 

e)  Isomach    Lines. 
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